#8/12/2015
#replication materials for Spirling "Democratization and Linguistic Complexity", JOP

rm(list=ls())

#Figure 9: proportion of all speeches are made by cabinet vs non cabinet
#Figure 10: probability a minister speaks following a non-minister for 'middle period'
# around claimed break
#(Figure 10 requires the 'effects' package)


############################
### Figure 9 ###############
############################


#load the data
#setwd("C:/Users/as9934/Dropbox/complexity/August2015/JOP_replication_data/")
load("bigframe.rdata")

sessions <- unique(big.frame$year.dummy)
speeches.cab <- c()
speeches.non <- c()
tot.speeches <- c()
mean.len.cab <- c()
mean.len.non <- c()


for(i in 1:length(sessions)){
  dat <- big.frame[big.frame$year.dummy==sessions[i],]
  tot.speeches <- c(tot.speeches, nrow(dat))
  speeches.cab <- c(speeches.cab,nrow(dat[dat$cabinet==1,]))
  speeches.non <- c(speeches.non,nrow(dat[dat$cabinet==0,]))
  
  mean.len.cab <- c(mean.len.cab, mean(dat$word.count[dat$cabinet==1]))
  mean.len.non <- c(mean.len.non, mean(dat$word.count[dat$cabinet==0]))
}

par(bg='cornsilk1')
par(mfrow=c(1,1))
plot(speeches.cab/tot.speeches,lwd=2, ylim=c(0,1), type="l", axes=F, xlab="", ylab="")
axis(1, at=1:length(sessions), labels=sessions)
axis(2)
lines(speeches.non/tot.speeches, lty=2)
legend("topright", lty=c(1,2), lwd=c(2,1), legend=c("cabinet","non"))
box()

############################
### Figure 10 ###############
############################

X11()
require(effects)

#function to regressions and plots
#(egregious use of <<- ... avert thine eyes)
make.pred <- function(sess="1865_1", pre = 0){
  big <<- big.frame[big.frame$year.dummy==sess,]
  curr <<- as.numeric( big$cabinet[-c(1)]) 
  prev <<- as.numeric( big$cabinet[-c(length(big$cabinet))] )
  mod <<- glm(curr ~ prev, family='binomial')
  out <<- as.data.frame(effect("prev", mod))[pre+1,]
  out
}

#do first one of each
preds0 <- make.pred(sess= "1852_1", pre=0) # this is Pr(cab now|bb)
preds1 <-  make.pred("1859_1", pre=1) # this is Pr(cab now|cab)


sessions <- c( "1852_2",  "1852_3",  "1852_4" , "1852_5" ,  "1857_2" , "1857_3"  , "1859_1","1859_2", "1859_3", "1859_4", "1859_5", "1859_6", 
               "1859_7", "1865_1", "1865_2", "1865_3", "1868_1", "1868_2", "1868_3", 
               "1868_4", "1868_5", "1874_1", "1874_2", "1874_3", "1874_4", "1874_5", 
               "1874_6", "1874_7")

for(i in 1: length(sessions)){
  preds0 <- rbind(preds0, make.pred(sessions[i], pre=0))
  preds1 <- rbind(preds1, make.pred(sessions[i], pre=1))
  
  
}






par(bg='cornsilk1')
plot(1:nrow(preds0), preds0$fit, ylim=c(min(preds0$lower),max(preds0$upper)), axes=F, 
     xlab="",ylab="", pch=16)
par(las=2)
axis(1,at=1:nrow(preds0), labels= c("1852_1",sessions))
axis(2)
for(i in 1:nrow(preds0)){
  arrows(i, preds0[i,4], i, preds0[i,5],code=3, angle=90, length=.05 )
}
box()

